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Abstract 

In the present study we investigate analytically the process of velocity and en- 
ergy relaxation in two-phase flows. We begin our exposition by considering the 
so-called six equations two-phase model [Ish75,Rov06]. This model assumes each 
phase to possess its own velocity and energy variables. Despite recent advances, 
the six equations model remains computationally expensive for many practical ap- 
plications. Moreover, its advection operator may be non-hyperbolic which poses 
additional theoretical difficulties to construct robust numerical schemes [GKCOl]. 
In order to simplify this system, we complete momentum and energy conservation 
equations by relaxation terms. When relaxation characteristic time tends to zero, 
velocities and energies are constrained to tend to common values for both phases. 
As a result, we obtain a simple two-phase model which was recently proposed for 
simulation of violent aerated flows [DDGIO]. The preservation of invariant regions 
and incompressible limit of the simplified model are also discussed. Finally, several 
numerical results are presented. 
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1 Introduction 



Currently, single phase flows modeling is essentially understood besides the 
problem of turbulence. In two-phase flows the situation is completely different. 
Nowadays there is no general consensus on the two-phase flow modeling. Tur- 
bulence modeling in two-phase flows is even more tricky [Llo05]. The two-phase 
models are based on space and time (ensemble) averaging of the local governing 
equations of each phase. Consequently, these models can only provide infor- 
mation on the average flow behaviour. The derivation of the proper model is 
far from being achieved [Ish75,DL79,DCL79,BD82,SW84,DW94,IH06]. One of 
the main difficulties lies in the determining the mass, momentum and energy 
transfers in the presence of steep gradients across the interfaces. Very often 
more or less accurate empirical correlations are used to describe such interface 
processes. In the same time, two-phase flows are very frequent in industry and 
in nature. Typical examples include water waves (especially during the wave 
breaking), petroleum and gas flow in pipes, nuclear reactors [SFW'^OS], etc. 

The most important industrial application of two-phase flows is the simulation 
of a wide spectrum of accident scenarios in Pressurized (or Boiling) Water Re- 
actors (PWR, BWR). Extensive experimental programmes for the PWRs are 
extremely expensive. Currently, new reactor design requires more and more 
numerical simulations of basic reactor features [SFW'''05]. The goal is to re- 
duce experimental programmes to a minimum. On the other hand, industrial 
codes (such as CATHARE [Bes90], RELAP5 [Ran85], THYC [ACOR95] or 
Neptune_CFD [BG05,GCMM09,GBBea07], for example) need to be supple- 
mented by a large number of empirical closures which prevents reliable code 
application outside of their validity domain. Also these codes use robust but 
higly diffusive numerical schemes which make it difficult the computation of 
strong variable gradients proper to many accident conditions. 

Two-phase models play also an important role in simulation of interfaces be- 
tween two immiscible compressible fluids. This problem is very challenging 
from numerical point of view. All existing methods to address this problem 
can be conventionally divided into four big groups: 

• Volume of fluid (VOF) [HN81,SZ99] 

• Level set method [OS88,OS94] 

• Lagrangian interface tracking [GGLTOO] 

• Diffuse interface method [Lar90,AMW98,Shy98] 

In the last approach we do not compute precisely the interface position. Due 
to numerical diffusion, the volume fraction takes intermediate values between 
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and 1 even if it is initialized by values in {0,1}. The resulting mixture of 
two fluids has to be properly modelled in the transition region. Two-phase 
models that we will discuss below provide a sound physical description of 
mixing zones. The diffuse interface approach is less accurate than mentioned 
above methods (unless some special sharpening techniques are employed) but 
it is much simpler to implement and can naturally handle arbitrary topological 
changes in the flow. 

Recently, a single velocity/ single energy two-phase model was proposed [Dut07] 
[DDG08a,DDG08b,DDG10] in ad-hoc way as a model for violent aerated flows. 
By analogy with the six equations model, this system was called the four equa- 
tions model since in one space dimension it consists of two equations of the 
mass conservation, one equation of the momentum and one of the energy con- 
servation. Two-phase models where both fluids share the same velocity were 
previously considered in [ACK02,Del05], for example. 

The main goal of the present paper is to show the connection between the 
six equations model [GKC01,GP05,Rov06,IH06] (two equations of mass, mo- 
mentum and energy conservation) and recently proposed four equations model 
[Dut07,DDG08a,DDG08b,DDG10] with single velocity and single energy. Phys- 
ically it is done by introducing relaxation terms into momentum and energy 
relaxation equations. These relaxation terms constrain velocities and energies 
to tend to a common value. Then, after taking the singular limit when the 
characteristic relaxation time goes to zero, one obtains a simplified model. 
Mathematically, it is achieved by employing a Chapman- Enskog type expan- 
sion [CC95]. This technique has been already successfully applied to the so- 
called Baer-Nunziato model [BGN86,BN86] in [MG05]. The present study is 
greatly inspired by their work. This type of reduction from the barotropic six 
equations model were recently done in [MDG09]. 

The present manuscript is organized as follows. We start our exposition by 
presenting the so-called six equations model in Section 2. Then, we complete 
this model by relaxation terms and derive in Section 3 a single velocity, single 
energy model. Preservation of invariant regions and incompressible limit of the 
resulting four equations model are studied in Sections 4 and 5 correspondingly. 
Then we present some numerical results in Section 6. Finally, this paper is 
ended by drawing out main conclusions of this study and some perspectives 
for future research (Section 7). 



2 Mathematical model 

Consider a fluid domain f2 C which is filled by two miscible fluids. All 
physical quantities related to the heavy and light fluids will be denoted by the 
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Figure 1. A quasi uniform air/water mixture after the wave breaking has occured. 

superscripts + and — correspondingly. For example, if we consider a plunging 
breaker with important mixing processes, then water velocity, density, viscos- 
ity, etc. will be denoted with the superscript +, while variables related to the 
air will have the superscript — . When the interfaces are so multiple that it 
is impossible any more to follow them (see Figure 1 for an illustration), the 
classical modelling procedure consists in applying a volume average operator 
[Ish75,Rov06]. At this level of description, two additional variables naturally 
appear. The so-called volume fractions t), x are defined as 

a^{x,t) := lim i^^, (1) 

the heavy fluid occupies volume dQ~^ C dQ and the light one the volume 
dQ~ C dfl (see Figure 2) such that 

\dn\ = \dn-^\ + \dn-\. (2) 

After taking the limit in relation (2), one readily finds 

a+{x,t) + a-{x,t) = 1, ^{x,t) eQ X [0,T]. (3) 

The volume fractions characterize the volume occupied by the correspond- 
ing phase per unit volume of the mixture. 

After applying the averaging procedure, we obtain the so-called six equations 
model [Ish75,SW84,GKC01,Rov06]: 

dtia^p^) + V ■ (a^p^u^) = 0, (4) 

dt{a^p^u^) + V ■ {a^p^u^ ® u^) + a^Vp = V ■ (a^r^) + a^p^g, (5) 
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Figure 2. An elementary fluid volume d^l occupied by two phases. 

dtia^p^E^) + V ■ (a^p^H^u^) + p dta^ = 

V ■ (a^T^u^) - V ■ (a^q^) + a^p^g ■ u^, (6) 

where t), u'^{x,t), E^{x,t) are densities, velocities and energies of each 
fluid respectively. The total energy is the sum of the internal and kinetic 
energies i?^ := + The specific total enthalpy is defined as 

:= + = E^ H — — , where is just the specific enthalpy. Finally, 

P 

g stands for the vector of the gravity acceleration. 

The symbol denotes the viscous stress tensor of each phase. If we assume 
both fluids to be Newtonian, the viscous stress tensor can be written as 

= X^trlD){u^)I + 2p^B{u^), trB{u^) = V ■ , (7) 

where I := (5ij)i<jj<3 is the identity tensor, 3(u) := ^(^u + *(Vw)^ is 

the deformation rate and A^, p^ are viscosity coefficients. For ideal gases, for 
example, these coefficients are related by Stokes relation + = 0. Recall 
that the derivation procedure presented in the next section does not require 
any particular form of the tensor r^. 

Finally, the heat flux in each fluid is denoted by (see Equation (6)). As 
with viscous stress tensor, we do not assume any particular form of the heat 
flux. However, in most cases the Fourier's law [Fou22] is adopted: 

:= -K^VT^, 

where is the thermal conductivity coefficient and is the thermodynamic 
temperature. 

Remark 1 When a material demonstrates flagrant anisotropic properties, it 
is advised to use the generalized Fourier's law [IH06]: 

q^ := ■ VT±, 

where is the conductivity tensor. The same remark applies to the viscous 
stress tensor as well. 
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In order to close the system (4) - (6), we have to provide two equations of state 
which relate the pressure p to other thermodynamic variables p = p^{p^, e^). 
We assume the equations of state to satisfy some general thermodynamic 
assumptions: 



.±/„± „±^ ^ n dp^{p^,e^ 



p±(p±,e±) >0, 



The symbol ^ , denotes the partial differentiation with respect to 



dp^ 



> 0, for p"^ > 



dpi 



the density p^ while the entropy is kept constant. Physically this quantity 
corresponds to the squared velocity of pressure waves in each fluid (see Lemma 
6). 

Remark 2 We reiterate that the model (4) - (6) is a one pressure model, 
i.e. two phases share the same pressure. Models with two pressures can also he 
derived [SW84,BN86,BGN86,RH88,MG05]. In this situation, an extra closure 
(algebraic or differential) must be added. In many cases this additional closure 
represents a relaxation process which will tend to equilibrate the pressures. 

Remark 3 The momentum balance Equation (5) can be rewritten differently: 

dt{a^p^u^) + V-(a^p^u^®u^ + a^pl)-pWa^ = V ■{a^T^) + a^p'^g. (9) 

These two forms are obviously equivalent for smooth solutions. However, it is 
not the case for discontinuous ones. We believe that this form is more relevant 
from physical point of view, since the flux of momentum involves the pres- 
sure. In our analytical investigations we consider smooth solutions, however 
for numerics we advise using the last form (9). 

Remark 4 While considering two-phase flows, it is useful to introduce several 
additional quantities which play an important role in the description of such 
flows. The mixture density p and mass fractions rrv^ are naturally defined as 

p(f,t) := a+p+ + a"/9" > 0, V(f , t) G x [0, T], (10) 

, a^p^ , _i_ _ ^ 

m := , consequently m + m =1. 

P 

The total density p is assumed to be strictly positive everywhere in the domain 
f2. Hence, the void creation is forbidden in our modeling paradigm. Important 
quantities p, rn^ will appear several times below. 

The six-equations two-phase model presented in this section contains some 
other simplifications. Namely, we neglect capillarity effects which could be 
taken into account in the form of the Korteweg term [Kor01,BDGG09]. More- 
over, terms modeling mass, momentum and energy exchange between the 
phases are neglected as well. In general, their form depends strongly on the 
flow regime under consideration and is a subject of current debates. 
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Remark 5 Even if the Jacohian matrix of the advection operator in (4) - (6) 
is not diagonalizable over R, the existence of global weak solutions was recently 
proven at least in barotropic case [BDGG09] if we supplement the model by 
viscous and surface tension effects [KorOl]. Previously, it was already shown 
that viscosity and thermal diffusivity are indeed sufficient for the evolution 
problem well-posedness in the linearized or in the nonlinear sense for small 
data [Ara80,Ram00]. 

Despite recent advances in numerical methods [GKC96,GKC01,Rov06] and 
several simplifications made above, the system (4) - (6) is still very complex 
to solve numerically for industrial scale problems. One of the main difficulties 
lies in the advection operator which is known to be generally non-hyperbolic. 
In the next section we will derive a simplified two-phase model which has fewer 
variables and possesses the requested hyperbolic structure for quite general 
equations of state. 



3 Relaxation process 



In this section we will supplement momentum and energy conservation Equa- 
tions (5), (6) by relaxation terms. Physically, they represent the friction or 
drag force between two phases. When time evolves, this mechanism will en- 
sure the convergence of two velocities and energies to common for both 
phases values u and E. Augmented Equations (5), (6) take the following form: 

dt{a^p^u^) + V ■ {a^p^u^ ®u^) + a^Vp = V ■ (a^r^) + a^p^g± F^, (11) 

dtia^p^E^) + V ■ (a^p^H^u^) + p dto^ = 

V ■ (a^T^u^) - V ■ (a^q^) + a^p^g ■ ± E^. (12) 

In the present study we choose friction terms and E^ in this form: 

^ K a-p-a+p+ ^+ ^_ K a-p-a+p+ ^ 

e a~^p+ + a;~p~ e a+p+ + a" p~ 

where u := mfu^ + m~u~ is the barycentric velocity, k = 0{1) is a dimen- 
sionless constant and £ is a small parameter which controls the magnitude of 
the relaxation term. Physically it represents the characteristic relaxation time. 
In the following we are going to take the singular limit as the relaxation pa- 
rameter e — )■ 0. This goal is achieved with a Chapman- Enskog type expansion 
[CC95]. 

This kind of computations has already been successfully carried out for the 
Baer-Nunziato model [BN86]. In this section we follow in great lines the work 
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of Guillard & Murrone [MG05]. However, the computation details are sub- 
stantially different. 



The first step consists in rewriting the governing Equations (4), (11), (12) 
in quasilinear form. To shorten the notation, we will also use the material 
derivative which is classically defined for any smooth scalar function (p{x,t) 
as 

d'^ip dip 



dt 



dt 



In computations below, we will need the following technical lemma: 

Lemma 6 Consider two compressible fluids with general equations of state 
P = P'^iP^i c"*^)- Then, partial derivatives of pressure p are given by: 



dp 



ds± 



dp 



dp 



1 



( p 



±\2 



X 



where s^, cf are the entropy and the sound velocity in the phase ± respectively. 



We also denote by x '■= 



and by 9^ := ^ 



PROOF. Obviously, one can write: 

dp 



dp 



dp^ 



dp 



ds 



± 



(13) 



Similarly, if one takes into account the definitions of and 9^: 

de^ = x^dp^ + d^dp. 
Now we will use the Gibbs relation [Ish75,Cal85,IH06]: 



de^ = T^ds^ + 



P 



±\2 



dp^ 



(14) 



(15) 



After expanding de^ and dp according to (14), (13) correspondingly, we get 
the following differential inequality: 



p 



rfp± + ( r± _ 



dp 



ds- 



ds^ = 0. 



Since this equality must be true for arbitrary de"^ and dp, we conclude the 
proof by requiring the coefficients to be equal to zero. 
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Theorem 7 Smooth solutions to Equations (4), (H), (12) satisfy the follow- 
ing system: 

a^p^T^—^ = a^T^ : Vu^ - V ■ (a^q^) ± {Ea - Ed ■ u^), (16) 
^±p±^^ + «±Vp = V ■ («±r±) + a^p^g± F^, (17) 

a^T^ : Vw^ - V ■ (a^q^) ± (^d - Fa ■ u^), (18) 
where is i/ie entropy, {cfY and 9^ are defined in Lemma 6. 



PROOF. The mass conservation Equation (4) is straightforwardly rewritten 
using the material derivative: 



d'^^a^p 
dt 



+ a^p^V-M^ = 0. (19) 



If we multiply the latter by il''^ and subtract the result from (11), we will 
obtain announced above Equation (17). Similarly, multiplying Equation (19) 
by E^ and subtracting it from (12) leads to the total energy equation in 
quasilinear form: 

d^E^ 

a p — h V • (a pu ) + paf = 

V ■ {a^T^u^) - V ■ {a^q^) + a^p'^g ■ ± E^. (20) 

The kinetic energy evolution equation is straightforwardly obtained after tak- 
ing the scalar product of (17) with «=*=: 



d 



'±rl|^±|2 



a^p^ — dt ^ ^^^^ ■ Vp = tt ^ ■ V ■ (a^r^) + a^p^g- u^±Ed-u^. (21) 



Recall that the total energy E"^ is constituted of internal and kinetic energies: 
E^ = e± + ll-u^p. Consequently, the internal energy evolution equation is 
derived by subtracting (21) from (20): 

a^p"^— VpV ■{a^u^)+paf = a^T^ : Vu^ -V ■(a^q^)±(Ed- Ed-u^). 

dt 

(22) 

In order to introduce the entropy variable into our considerations, we will 
make use of the Gibbs relation. After multiplying (15) by a^p^ and passing 
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to the material derivative, one gets: 



dt 



dt 



a 



'-- 'dF' 



Substituting (22) into the last relation and taking into account (19), leads to 
the requested entropy Equation (16). 

From lemma 6 we directly obtain the following equality: 

After substituting (16) and replacing (^jj^ — X^) by 6^{cf)'^, we obtain (18) 
and the proof is completed. 



Computations that we will perform below will be clearer if the governing 
Equations (16) - (18) are recast in the vectorial form: 

A{V,)^ + B{V,)VV, = V ■ T{V,) + S(y,) + (23) 

where we introduced several notations. The vector represents six unknown 
physical variables := *(s^, s~, -u"*", m a"*") and ^ is the componentwise 
partial time derivative. Symbol VVJ, denotes this vector: 

:= *(Vs+,Vs-,(-V)ii+,(-V)n-,Vp,Va+). 
Matrices A(Ve) and B(K) are defined as 



MV.) := 



(4 

















\ 
























a 




-I 



















p" 


-I 



















a+ 






u 













-Pt 
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10 



B(V;) := 









































n 


n 


rv^ it ~^ 

U-t jJ (Jj 


n 


a+l 


n 















a^l 




















pju + 




V 












-PqU- 


) 



where several symbols were introduced to shorten the notation: 

p± := (p±)2(c±)26»±, := a^p^T^, := a^p^O^. 

In these matrix notations the size of zero entries must be chosen to make the 
multiplication operation possible. 

On the right hand side of (23), the work of viscous forces is denoted by symbol 
V ■ TiVe) which is defined as 

V ■ T(V;) := *(a+T+ : Vn+ - V ■ (a+g+), a"r" : ViT" - V ■ 
V-(a+T+), V-(a"T"),Q;+r+ : Vm + -V- (a+g"*"), a'r" : Vu'-V-{a'q')). 

The source term SiVe) := *(0, 0, a^p^^, 0, 0) incorporates the gravity 

force and R(VI;) contains the relaxation terms: 



R(K) := K- 



a^p'^a p 
a+p+ + a~ p~ 



1, — 1, -u — M — (m — M )^ X — M ). 



Since we expect the limit -^V to be finite as e — > 0, necessary the limiting 
vector V lies in the hypersurface R(l^) = 0. In terms of physical variables, it 
implies u = = u~ . Consequently, we find our solution in the form of the 
following Chapman-Enskog type expansion [CC95]: 

V, = V + eW + 0{e^). 



After substituting this expansion into (23) and taking into account the fact 
that R(V^) = 0, at the leading order in e one obtains: 

dV 

A(V^)^ + B{V)VV = V ■ T{V) + S{V) + R'{V)W, (24) 
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where 



R'{V) 



K- 



+ a~p' 



^0 








0^ 

















I 


-10 





-I 


I 


























Oy 



Henceforth, we make a technical assumption of the presence of both phases in 
any point x G i7 of the flow domain. Mathematically it means that < < 1. 
Since + = 1, obviously the same inequality holds for a~ . Otherwise, 
the relaxation process physically does not make sense and we will have some 
mathematical technical difficulties. 

Under the aforementioned assumption, the matrix A{V) is invertible. Hence, 
we can multiply on the left both sides of (24) by PA~^(V^) where the projection 
matrix P is to be specified below: 



dV 

P— + PA-^(y)B(y)Vy = PA-\V)V-T{V) + PR'{V)W+PA-\V)S{V), 



where R'{V) := A ^{V)T{.'{V) and has the following components 



(25) 



K'iV) 



K- 



a^p^a p 














1 



a'^p'' 



a p 










-10 



a^p 




10 






The vector of physical variables V has six components (in ID): 

V = ''{s^ , s~ ,u ,u ,p, a'^) 

and only five are different. In order to remove the redundant information, 
we will introduce the new vector U defined as U := ^{s~^ , s~ ,u ,p, a^). The 
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Jacobian matrix of this transformation can be easily computed: 

/l o\ 



J := 



dV 
dU 



10 

10 

10 

1 

yo ly 



In new variables Equation (25) becomes: 
dU 

PJ—+PA-\U)B{U)3VU = PA-\U)V-T{U)+PR'{U)W+PA-\U)S{U). 

(26) 

Now we can formulate two conditions to construct the matrix P. First of 
all, the vector W is unknown and we need to remove it from Equation (26). 
Hence, we require PR'(\^) = 0. Then, we would like the governing equations 
to be explicitly resolved with respect to time derivatives. It gives us the second 
condition PJ = I. The existence and effective construction of the matrix P 
satisfying two aforementioned conditions 

PR'{V) = 0, PJ = I, 

are discussed below. Presented in this section results follow in great lines 
[MG05]. 

We will consider a slightly more general setting. Let be vector G and its 
reduced counterpart U G M""'^, A; < n. In such geometry, R'(l^) G Mat„,„(R), 



J e Mat 



n,n—k 



and, consequently, P G Mat 



„_A;,„(]R). Here, the notation 
Matm,n(l^) denotes the set of m x n matrices with coefficients in R. We have 
to say also that from algebraic point of view, matrices R'(^) and R'(^) are 
completely equivalent. Thus, for the sake of simplicity, the following proposi- 
tions will be formulated in terms of the matrix R'(V). 

Lemma 8 The columns of the Jacobian matrix J form a basis of ker(ll'{V)^ . 



PROOF. If we differentiate the relation R(V) = with respect to U, we 
will get the identity 'R'{V)3 = 0. It implies that range(^J^ C ker(^R'(y)^. 

By direct computation one verifies that dim range (R'(y)) = k. From the 

well-known identity range (R'iV)) © ker(R'(V)) = M", one concludes that 
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dimker(^R'(l^)j = n — k. But in the same time, the rank of J is equal to n — k 
as well. It proves the result. 

Theorem 9 We suppose that for all V, range (R'(y)) n ker(R'(y)) = {0} 
then there exists a matrix P G Mat„_A;^„(]R) such that PR'(y) = and PJ = 

Ira— fc ■ 



PROOF. Assumption range (R'(y)) nker(R'(y)) = {0} implies that 
range(R'(y)) ©ker(R'(F)^ 



From Lemma 8 it follows that range (j) = ker(R'(y)). Thus, the space M" 
can be also represented direct sum range (R'iV)) ©range(j). We will 
define P to be the projection on ker(^R'(y)^ = range(^J^. Since obviously 
R'(y) G range (^R'(V^)^ and J G range(^J^, we have two required identities: 
PJ = In-k and PR'(y) = 0. 



Now, in order to compute effectively the projection matrix P, we will con- 
struct an auxiliary matrix D(y) = [J^, . . . , J""'^, , . . . , I''], where J* is the 
column i of the matrix J and {/^, . . . , /^} are vectors which form a basis of 
range(^R'(y)^ . We remark that PD(V^) = [l„_fc,0]. Lemma 8 implies that 
the matrix D(l/) is invertible. Thus, the projection P can be computed by 
inverting Y){V): 

P=[In-k.,0]-B~\V). (27) 



Let us apply this general framework to our model (25), where n = 6 and k = 1. 
The matrix D(y) and its inverse D~^{V) take this form: 



'^l ^ 
10 
1 m-I 
1 -m+I 
10 

yo 1 J 

where are mass fractions defined in Remark 4 

















o\ 





1 






















m~I 




















1 




















1 







I 


-I 
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Now, the projection matrix P can be immediately computed by (27): 

^1 0^ 
10 
m+I m~I 
10 
1 



Finally, after computing all matrix products PA ^(t/)B(t/)J, PA ^(^7)V ■ 
T{U), PA"^(f/)S(f/) present in Equation (26), we obtain the following model: 



du 



~dt 



a^T^ : Vw=^ - V ■ (a*g 



p— + Vp = pg + V -T, 



(28) 
(29) 



U 



P (c; 



9+p^ 



■(^a+r+ : - V ■ (ct+f 



+ (a-T- : - V ■ (a-q-)) , (30) 



e-p- 



n( 



da^ 
~dt 



+ a~^a SV ■ u 



a 



e+p^ 



a^T^ : Vu — V ■ {a^q' 



a-T' -.Vu -V ■ {a'q')), (31) 



where p = a"'"p+ + a p is the mixture density (10), physically is the sound 
velocity in the mixture: 

^pc^=p+(c+)V(c;)^ (32) 

To shorten the notation, we also introduced a quantity 6 defined by 

U6:=p^ctr~p-{c:r (33) 
and n := a'p+{c^y + a+p"(c7)^. 



The viscous stress tensor of the mixture is defined as r := a~^T~^ + a r .If 
both fluids are assumed to be Newtonian (7), we can express it in closed form: 



T := XtTB{u)I + 2pJ]){u), A:=a+A+ + a A , /i := a+^+ + a p . 



15 



Similarly, the thermal flux q in the mixture is given by this formula: 

q := a^q^ + a^q^ . 

For practical computations it is better to rewrite Equations (28) - (31) in 
conservative form [Lax73,GR96] which is valid for discontinuous solutions as 
well. After repeating in inverse sense computations from Theorem (7) one 
obtains the following system: 

dt{a^p^) + V ■ {a^p^u) = 0, (34) 
dtipu) + V ■ (pu +pl) = pg + V ■ T, (35) 
dt{pE) + V ■ (pHu) = pu -g + V ■ (tu) - V ■ g . (36) 

If we neglect viscous stress r and thermal flux q, we will obtain the so- 
called four equations model recently proposed by F. Dias, D. Dutykh and 
J.-M. Ghidagha [Dut07,DDG08b,DDG10,DDG08a] for simulation of violent 
aerated flows. Formal computations presented in this study can be considered 
as a step towards justiflcation of the four equations model. It can be also shown 
[DDGIO], that the advection operator of Equations (34) - (36) is hyperbolic 
for quite general equations of state (8). Hence, derived here four equations 
model is very attractive from modeling and computational point of view. 



4 Invariant regions 



Consider only the hyperbolic part of the system (34) - (36) together with 
source terms due to gravity: 

dt{a^p^) + V ■{a^p^u) = 0, (37) 

dtipu ) + V ■ (pw ® M + pi) = pg, (38) 

dt{pE) + V ■{pHu) = pu -g. (39) 



In this section we would like to study the preservation of invariant regions 
under the dynamics of the system (37) - (39). This property is very important 
in practice for stability of computations. Namely, from physical sense and from 
deflnition (1) of volume fractions it follows that < < 1, V(x, t) G 
Q X [0,T]. The question we address here is whether the region < a± < 1 
remains invariant under the system (37) - (39) dynamics. This property can be 
checked only for a"*", for example, since + a~ = 1. If the answer is negative, 
the model can be hardly applied to any practical situation. Positive result was 
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already proved for a barotropic model with two velocities in [BDGG09]. In the 
same work of Bresch et al. (2009) it was also shown that two-phase mixtures 
of compressible/incompressible fluids do not preserve the invariant regions. 

The notion of invariant regions is due to Chueh, Conley and Smoller [CCS77]. 
Independently it was also introduced by Weinberger [Wei75] for scalar parabolic 
equations. 

In order to handle this problem, it is more convenient to rewrite Equations 
(37) - (39) in terms of the so-called physical variables W =* {u ,p,a~^ , s). 
The derivation of the equation for the velocity u evolution is straightforward. 
Equations for thermodynamical variables p, and s are trickier to obtain. 
However, it was already done in [DDG08a,DDG10] and partially above (see 
Equations (28) - (31)): 

Lemma 10 Smooth solutions to Equations (37) - (39) satisfy the following 
quasilinear system: 



dtu + u - Vu + — = g, (40) 
P 

dtp + u ■Vp + pclV -u = 0, (41) 

dta~^ + u ■ Va+ + a+a~5V ■ m = 0, (42) 

9tS + M-Vs = 0, (43) 

where pc^ and 5 are defined in (32) and (33) respectively. 



PROOF. For the proof see, for example. Bias et al. (2009) [DDGIO] where 
the authors worked with an auxiliary variable a := — a~ . Taking into 
account the relation + a~ = 1, one can express the volume fractions in 

terms of a: = — - — . Consequently, a'^a^ = — - — . 



Corollary 11 System (40) - (43) which governs the evolution of physical 
variables can be recast under the following matricial form: 

dW dW 

^ + M(W^)^ = S(1^)> (44) 
where W = {u ,p, a"*", s), source terms S{W) = ^{g, 0, 0, 0) and ifn = (ni, n2, n^) 
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is a normal direction, Un = u ■ n is a normal velocity, then 

( 



m{W)n = Y,Mi{W) 



rii 



i=l 
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Now we can state the main result: 

Theorem 12 For smooth solutions to Equations (37) - (39), the region < 
<1 remains invariant under the system dynamics. 



PROOF. The proof is based on the theory of J. SmoUer [Smo94]. We know 
that the value a" = corresponds to W5 = a+ = 1. Obviously, the case 

= Q can be treated similarly with the same conclusions. The matrices 
MiiW) are smooth functions of W in the neighborhood of W which satisfies 

G [0, 1]. We would like to check whether the boundary = 1 is invariant 
under the dynamics of (44). According to the theory of invariant regions by 
Chueh, Conway and Smoller ([CCS77]), this will be the case if and only if 
d{W^ — l) = dW^ = (0, 0, 0, 0, 1, 0) is a left eigenvector of matrices Mi{u ,p, 1, s) 
for all admissible values of u , p and s. By straightforward computation one 
can easily check that 

(0, 0, 0, 0, 1, 0) ■ Mi{u,p, 1, s) = Ui{0, 0, 0, 0, 1, 0). 

The theorem is shown now. 



Remark 13 The present result shows also that if negative values of volume 
fractions are reported in simulations, they are due to some numerical in- 
stabilities and have nothing to do with mathematical properties of the four 
equations model (37) - (39). 



5 Incompressible limit 



In practice, there are many situations when governing equations can be fur- 
ther simplified by filtering out acoustic effects. For various combustion models 
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it was done by A. Majda and his collaborators [Maj82,KM82,MS85]. More re- 
cently these techniques were applied to some two-phase models [KPTV00,Del05], 
[DelOT]. 

In the present section we derive the incompressible limit of the four equations 
model (34) - (36). For the sake of simplicity we will neglect thermal fluxes 
which are not important for the exposition below and do not directly affect 
acoustic waves propagation. Hence, in this section we consider the following 
set of equations in nonconservative form for convenience: 



dt{a^p^) + V ■ {a^p^u) = 0, (45) 
pdtu + p{u ■ V)u + Vp = pg + V ■ T, (46) 
pdte + pu ■ Ve + pV ■ u = t : Vu , (47) 

where we expressed the energy balance in terms of the internal energy e. 

In order to estimate the relative importance of various terms, we introduce 
dimensionless variables. In this section we use the special low Mach number 
scaling which reveals acoustic effects magnitude. The characteristic length, 
time, and velocity scales are denoted by £, and Uq respectively. For example, 
I can be chosen as a linear dimension of the fluid domain f2. The density, 
viscosity and sound velocity scales are chosen to be those of the heavy fluid, 
i.e. po 5 ^"^^ correspondingly. Since we are interested in acoustic effects, 
the natural pressure scale is given by pt{cQsY- If 'we summarize these remarks, 
dependent and independent dimensionless variables (denoted with primes) are 
defined as: 



- t' ■=- io^y ■= — iu^y ■= 

<^ ^0 Po Po ^0 



■ - !L ' = P 

Since the volume fraction is dimensionless by definition (1) we keep this vari- 
able unchanged. 

After dropping the tildes, nondimensional system (45) - (47) of equation be- 
comes: 

St dt{a^p^) + V ■ (a^p^u) = 0, (48) 
St pdtu + p{u ■ V)u + ^ Vp = ^pg + ^ V ■ T, (49) 

St pdtC + pu ■'Ve+ ■ u = -^r : Vu , (50) 

Ma Re 

where several scaling parameters appear: 
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i 

• Strouhal number St := . In this study we will assume the Strouhal 

Uoto 

g 

number to be equal to one St = 1, i.e. to = -rr- 

U 

• Mach number Ma := which measures the relative importance of the 
flow speed and the sound speed in the medium. 

• Proude number Fr := ° compares inertia and gravity forces. 

U i 

• Reynolds number Re := — qr gives the measure of the ratio of inertial to 
viscous forces. 

In this section we will consider the asymptotic limit as Ma — )■ 0. Consequently, 
all physical variables a"^, p^,p,u and e are expanded in formal series in powers 
of the Mach number: 

ip = (fo + Ma(pi + Ma'^ip2 + . . . , (p E {a"^, p'^,p,u ,e}. (51) 

Formal expansion (51) is then substituted into the system (48) - (50). Our 
goal is to derive a system which governs the evolution of physical variables 
at the lowest order in the Mach number (e.g. -Uq, a^, etc.). Equation (49) at 
orders Ma~^ and Ma~^ leads: 

Vpo = Vpi = 0. 

In other words, po = Po{t) and pi = Pi{t) are only functions of time. Internal 
energy balance equation (50) gives us the incompressibility constraint V -Uq = 
at the leading order Ma~^ and, correspondingly, V-m i = at the order Ma~^. 

Relation (3) after asymptotic expansion (51) becomes 

where Sok is the Kronecker delta symbol equal to 1 if A; = and otherwise. 
Thus, at the leading order in Mach number we keep usual relations: 

"o' + «o=l' Po = OoPo + «oPo • 
Equations (48), (49) at the order Ma° become: 

d,{a^p^) + V ■ (a^p^uo) = 0, (52) 



podtUo + po{uo ■ V)uo + Vtt = :^Pog + —V ■ tq, (53) 



where by tt we denote the pressure oscillations p2 at the order Ma^. 
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Finally, some information on the behaviour of can be deduced from the 
Gibbs relation (15) which takes the following dimensionless form: 



T^ds^ = de^ ^?r^dp^- 

Ma2(p±)2 

After dividing by dt and expanding the Gibbs relation as in (51) we get these 
results: 

^ = 0, ^ = 0. (54) 
dt ' dt ^ ' 

In particular, it means that if in each phase the density p^^ was initially 

constant ^ , it remains so under the system dynamics. Moreover, taking into 

account the first relation in (54) and the flow incompressibility constraint 

V ■ Mo = 0, one can easily deduce from the mass conservation Equations (52) 

that volume fractions are also simply transported by the flow: 

dtUQ + uo ■ Va^ = 0. 



If we summarize all the developments made above and turn back to physi- 
cal variables, we will get the following system of equations which governs an 
incompressible two-phase flow: 



dta^ + u ■ Va^ = 0, (55) 
V ■ u = 0, (56) 
pdtu + p{u ■V)u +^7^ = pg + V ■ T, (57) 

where we dropped the index to simplify the notation. The mixture density 
p and the viscous stress tensor r are defined as above. This system should be 
completed by appropriate boundary and initial conditions. 

Remark 14 Two- fluid incompressible models, such as the system just derived 
above (55) - (57), are often used in many practical situations such as wave 
breaking [CKZL99] and others [SZ99]. For powder-snow avalanches applica- 
tions, the volume fraction Equation (55) is completed by a diffusive term ac- 
cording to the Pick's law: 

dta^ + u ■ Va^ = V ■ (z//Va^). 

This extra term allows to take into account the mixing between two fluids due 
to turbulence effects, for example. For more details on this extension we refer 
to [JR93,ESH04,DARB09J. 

^ This assumption is often made in applications. For example, pure water and pure 
air are assumed to have constant densities under normal conditions. However, the 
mixture density can undergo strong variations (e.g. after the wave breaking). 
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6 Numerical results 



In the present section we will consider only the advective part of the system 
(34) - (36), i.e. we neglect viscous and thermal fluxes: 

dtia^p^) + V ■ (a^p^u) = 0, (58) 
dt{pu) + V ■ {pu +pT) = pg, (59) 
dtipE) + V ■ (pHu) = pu ■ g. (60) 

Actually, the system we consider here is stiffer than the original Equations 
(34) - (36) since it does not contain any diffusive effects. The resulting equa- 
tions have been shown to be hyperbolic for any reasonable equation of state 
[DDG10,DDG08a]. For illustrative purposes we assume here that both fluids 
are governed by stiffened gas type laws: 

where 7^, vr^ and are some constants determined by physical properties 
of pure fluids under consideration. Also we assume that two fluids are in the 
thermodynamic equilibrium: 

p = = , T = = T~ . 

We refer to [Dut07,DDG10,DDG08a] for more details on the construction of 
this equation of state and the discussion of some properties of such two-phase 
mixtures. 



6.1 Numerical schemes 



For the numerical study we choose the finite volumes method [Kro97,BO04] 
since it is the method of choice for the systems of conservation laws due to its 
excellent local conservative properties. More precisely, we use the cell-centered 
approach [BJ89,Bar94] which is more natural in our opinion. For simplicity we 
assume that the system of equations is solved in M^. However the extension 
to 3D for cartesian meshes is straightforward. We briefiy describe below the 
discretization procedure adopted in this study. References are also provided if 
more details of the discretization method are needed. 
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6.1.1 Space discretization 

System (58) - (60) can be written as 

^ + V-^(w) = S(w), (61) 

where 

w = (?yi)f=i := (a+p+, a^p", pwi, pM2, pE) , 
and, for every n = (ni,n2) G M^, 

J^(w) ■ n = {a^p'^u ■ n, a~p~u ■ n, pu ■ nui + prii, 

pu ■ nu2 + pn2, pHu ■ n) , (62) 

S(w) = (0, 0, pgi, pg2, P9 ■ u) ■ 
Then, the Jacobian matrix A(w) ■ n is defined by 

. / N . / N ^ d(J^(w) ■ n) , ^ 

A„ w := A w • n = -^-^ ^ . 63 

In order to compute A(w) ■ n, one writes Equation (62) for J^{w) ■ n in terms 
of w and p: 

^, . ^ ( w^rii + Win2 wsUi + Win2 wsrii + W4n2 
J^{w)-n=[wi ■ , W2 ■ , W3 ■ \-pni, 

\ W1 + W2 Wi + W2 Wi + W2 

Wsni + W4n2 , .Wsni + W4n2\ 
W4 ■ \-pn2, {W5+P) ■ . 

Wi + W2 W1 + W2/ 

The Jacobian matrix (63) then has the following expression: 
A(w) ■ n = 

Un^^ -Un^ ^rii ^n2 

II p "■ p p ^ p ^ 

-Un^ ^n2 

-U2Un + -U2U^ + ^Tls ti2^1 + ^^^2 + ^^2^12 + ^^2 ^^2 

^nn{it,-H) u^{-t;-H) u^^^+Hn, u^^^+ Hn2 u^{l + ik 
where Un = u ■ n. 

The computational domain C is decomposed into a set of control volumes 
T such that Vt = Uk^tK. We integrate Equation (61) on K: 

^fwdn+y I J^{w)-nKLda= [ S(w)(i^], (64) 
at JK L£M{K) "'^n^ •'^ 
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where ukl denotes the unit normal vector on K n L pointing into L and 
Af{K) = {LeT : area(K n L) ^ 0} . Then, setting 



vol A JK 



Y0\{K) Jk 

we approximate (64) by 

dwK , ^ area(Lnir) ^ 

where the numerical flux 

^{wk, wl; Ukl) ~ ,1 ^ [ -7^(w) • Hkl da 

area(L fl K) Jkhl 

is explicitly computed by the FVCF formula of Ghidaglia et al. [GKCOl]: 
J'(wi^) ■ n + J^(wl) ■ n 



- sign(A„(/i(wii-, Wi))) . (65) 



Here the Jacobian matrix A„(/i) is defined in (63), /i(wi^, w^) is an arbitrary 
mean between wk and and sign(M ) is the matrix whose eigenvectors are 
those of M but whose eigenvalues are the signs of that of M. In this section we 
did not deal with boundary conditions. We refer to [GP05] for more details. 



6.1.2 Higher order extension 

In the previous section we described the first order scheme which might be 
too diffusive for most practical applications. That is why we present here 
a higher-order extension which is a variant of MUSCL^ limiting technique 
[Kol75,Kol72,vL79,vL06]. This numerical method ensures stability and non- 
oscillatory behaviour of numerical solutions. To describe this scheme we switch 
to Cartesian notation since a bigger stencil is needed for the gradient recon- 
struction procedure. In this notation Wj denotes the average of conservative 
variables in the cell i and w^'^ denotes respectively reconstructed left and 

right states at cell faces. According to the method adopted in the current 
study [Bar94,BO04,TMD07], cell faces values are computed as: 

wf+i =w, + U{l- i^Mrf ){w, - w,_i) + (1 + (w,+i - w,)) , (66) 



^ Acronym MUSCL stands for Monotone Upstream-centered Scheme for Conserva- 
tion Laws [vL79] 
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1 = Wi+i-^(^(l-/t)^(rf)(wi+2-Wi+i) + (l+«;)(/j(^)(wi+i-w,)^, (67) 
where k G [— 1, 1) is a free parameter and 



1+ 



Wi+i -Wi ^ W,;+i - Wi 



Wi - Wj_i Wi+2 - Wj_i 

Then, reconstructed values are used to compute the numerical flux (65) 
of the FVCF scheme. 

The function ip{r) is called the limiter function and is incorporated to obtain 
non-oscillatory resolution of discontinuities and steep gradients [BB73,Swe84]. 
We tested in practice two following limiter functions: 

• M3 limiter proposed in [ZD98,TMD07]: 

2Nr \ / 2r \^ 
l + r2/V 1 + r 



^M3(r) = l-(l + ^^ 1-— , N = 2, {61 



• MinMod limiter with compression parameter [Bar94,BO04]: 

'/'MM (r) = max{0,min{r,^}}, /3g(1,- . (69) 

In the smooth unlimited form {ip{r) = 1), the truncation error e for a scalar 
conservation law with the smooth flux f{u) is given by [Bar94,BO04]: 



Consequently, for k = |, the scheme provides theoretically in smooth regions 
an overall spatial discretization with 0{Ax^) error. In computations presented 
below, we use the optimal choice of this parameter together with M3 limiter 
(68). 

Remark 15 If the limiter function is symmetric, i.e. 

the reconstruction formulas (66), (67) are independent of the parameter k and 
degenerate to the classical MUSCL2 scheme. It is straightforward to show that 
limiter functions (68) and (69) are not symmetric. 

6.1.3 Time discretization 

In the previous section we briefly described the spatial discretization procedure 
with some finite- volume scheme. It is a common practice in solving time- 
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dependent PDEs to first discretize the spatial variables. This approach is called 
a method of lines: 



dw FV 

— + V ■ J'(w) = S(w) =^ wt = C{w), 

where C{w) is a space discretization operator. In order to obtain a fully dis- 
crete scheme, we have to discretize the time evolution operator. In the present 
work we decided to retain the so-called Strong Stability-Preserving (SSP) 
time discretization methods described in [Shu88,GST01,SR02]. Historically 
these methods were initially called Total Variation Diminishing (TVD) time 
discretizations. However, this term is not completely correct. In computa- 
tions presented below we use the following third order four-stage SSP-RK(3,4) 
scheme with CFL = 2: 



wW = w(") + ^At£(w(")), 

w(2)=w« + ^At£(w«), 

w(3) = + iw(2) + iAt£(w(")), 

3 3 6 ^ ^' 

w("+i) = w(3) + iAt£(w(=^)), 
where w^"^ = w(-,t„) and w*^""*"-^) = w(-,t„_|_i). 

The time step is chosen adaptively to satisfy the following stability condition 
[CFL67]: 

^ CFL ■ AXmin 

where Axmin is the minimal mesh spacing, u is the velocity and Cg is the sound 
speed in the mixture (32). 

Remark 16 Presented above time step restriction is used together with the 
first order finite volumes scheme. The application of the MUSCL3 scheme 
requires the following modification of the CFL condition: 

At < 



2 K maX-fjlt -|- Cslmax) |^ C^lmaa;} 

6.2 Two- fluid Sod shock tube problem 



We present here a two-fluid generalization of the classical Sod shock tube 
problem [Sod78]. The sketch of the initial condition is given on Figure 3. The 
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£ = 1 





a , p , p 




u = 


u = 







L = 2 

Figure 3. Sketch of the initial condition for the Sod shock tube test case. 





Velocity, u 


Density, 


Volume fraction, 


Pressure, p 


Heavy fluid, + 


0.0 


1.0 


0.98 


1.0 


Light fluid, — 


0.0 


0.125 


0.02 


0.1 



Table 1 

Initial condition parameters for the two-fluid Sod shock tube problem. 





7^ 


TT^, Pa 


' kg-K 


Heavy fluid, + 


2.6 


0.0 


661.0 


Light fluid, — 


1.4 


0.0 


661.0 



Table 2 

Equation of state parameters for the two-fluid Sod shock tube problem. 



Two-fluid Sod problem. Density plot Two-fluid Sod problem. Densify plof 




0.5 1 1.5 2 0.5 1 1.5 2 



(a) FVCF l^t order (b) MUSCL 3 

Figure 4. Convergence of the solution with /i-refinement 

gravity force is not taken into account in this test case, i.e. g = 0. Values of 
all parameters used in this computation are given in Tables 1 and 2. 

The simulation was stopped at time T = 0.4 s and the computation results 
are presented on Figures 4 and 5. On all Figures we show the total density p 
plots and perform a comparison between FVCF and MUSCL3 schemes. Figure 
4 shows the behaviour of the solution when the mesh is refined. Left image 
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Two-fluid Sod problem. Density plot: N = 200 Two-fluid Sod problem. Density plot: N = 400 




0.5 1 1.5 2 0.5 1 1.5 2 



(a) N = 200 (b) N = 400 

Figure 5. Convergence of the solution with p-refinement 

shows the results of the first order FVCF scheme, while the right one refers to 
the MUSCL3 extension. The improvement of the solution is obvious in both 
cases. In particular, we would like to underline the excellent shock resolution 
by the higher order scheme (Figure 4 (b)). 

We perform also the comparison between two schemes (see Figure 5) for the 
same mesh size. As it was expected, higher order scheme greatly improves 
the sharp resolution of discontinuities for nearly comparable CPU time. Thus, 
there is an obvious interest in using MUSCL type schemes for practical appli- 
cations. 

This test case clearly shows the convergence of the numerical solution when the 
mesh is refined (/i-convergence), but also with respect to the scheme order {p- 
convergence). Consequently, it validates our numerical code. In all subsequent 
computations we use the MUSCL3 scheme with the optimal choice of the 
parameter k = |. 



6.3 Water drop test case 



The sketch of this numerical experiment is given on Figure 6. The values of 
parameters are given in Table 3. Initially the velocity field is taken to be zero 
and the pressure field is uniform in all domain po = 10 Pa. The gravity force 
is taken to he g = 10 m/s^ and the computational domain is discretized with 
a uniform grid of 100 x 100 control volumes. Results of this simulation are 
presented in Figures 7 - 10. A similar test case has already been considered 
in [Dut07,DDG10,DDG08a] and results are in good qualitative agreement. 
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{xo,yo) 



a ,p 



Figure 6. Sketch of the initial condition and computational domain geometry for the 
water drop test case. In our simulation we took the following values of geometrical 
parameters: i = 1.0 (the domain is square), ro = 0.15, xq = 0.5 and yo = 0.7. The 
gravity acceleration is taken to be 5 = 10 m/s^. 





7^ 


TT^, Pa 


'^v > kg-K 


± kg 


a+ 


Heavy fluid, + 


1.6 


0.0 


1.0 


5.0 


0.99 


Light fluid, — 


1.4 


0.0 


1.0 


1.0 


0.01 



Table 3 

Equation of state parameters for the water drop test case. 




(a) t = 0.1 s (b) t = 0.2 s 



Figure 7. Water drop test case: initial acceleration and deformation stage. 
7 Perspectives and conclusions 

In this study we considered several two-fluid models. We began the exposi- 
tion by the so-called six equations model (4) - (6). Despite recent progress 
[ACOR95,GKC96,TKP99,GKC01,BG05,SFW+05,BQ06,Rov06,NKDVLG08], 
[GCMM09], this system still represents some major difficulties for the nu- 
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(a) t = 0.3 s (b) t = 0.4 s 

Figure 8. Water drop test case: falling and further deformation of the drop. 




(a) t = 0.5 s (b) t = 0.6 s 

Figure 9. Water drop hitting the tank bottom. 




Density 




(a) t = 0.7 s (b) t = 0.8 s 

Figure 10. Water drop flow on the bottom. 

merical solution. Namely, the advection operator may be non-hyperbolic and 
contains non- conservative terms to be defined in some sense for discontinuous 
solutions. That is why, the six equations model was simplified through the ve- 
locity and energy relaxation process (see Section 3). In this way, we formally 
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derived the so-called four equations model (34) - (36) recently proposed as a 
model for violent aerated flows [Dut07,DDG08b,DDG10,DDG08a]. Thus, the 
present work can be considered as an attempt towards further comprehen- 
sion and at least formal justification of single velocity, single energy two-phase 
models. The resulting system (34) - (36) is hyperbolic for any reasonable equa- 
tion of state [DDGIO] and possesses several nice properties. In particular, in 
Section 4 we show that invariant regions G [0, 1] are preserved under the 
system dynamics. This property is necessary for the well-posedness of the sys- 
tem. In Section 5 we also formally derived the incompressible limit as the Mach 
number tends to zero. As a result, we recover usual two- fluid incompressible 
Navier-Stokes equations [PZ99,SZ99,CKZL99] if both fluids are assumed to 
be Newtonian. Finally, several numerical results are presented in Section 6. 
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